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The material derivative idea of continuum mechanics (Ref. 1) and the 
adjoint variable method of design sensitivity analysis are used to obtain a 
computable expression for the effect of shape variations on measures of 
structural performance of three-dimensional elastic solids (Ref. 2). 

Consider the three-dimensional elastic solid shown in Figure 1, with the 

° i 7 2 3 T 

shape of the domain Q as a dgsigijt variable. In Figure 1, z = [z , z , z ] Ls 
the displacement field and T , T , and i are clamped > traction free, and 
loaded boundaries, respectively. 

Using the principle of virtual work, the variational equilibrium equation 
for the elastic solid can be obtained (Ref. 3), where a^J(z) and £*j(z) are 
the stress tensor due to a displacement z and the strain tensor due to a 
kinematically admissible virtual displacement z, respectively, f = [f , f , r ] x 
is the body force, T = [T , T , T*]™ is the traction force, and Z is the space 
of kinematically admissible virtual displacement. When the Galerkin method is 
applied to the variational equilibrium equation for approximate solution, an 
approximate finite-element equation is obtained. 


THREE DIMENSIONAL ELASTIC SOLID 



• Principle of Virtual Work: 

a (z,z) = JJ/J l o ij (z)E i j(z)]dn 
i>j=l 

- J/U l // 2 [ l T^Jdr- yz), 

“ i=l ? i-1 

for all z £ Z 

• FEM Equation is an approximate equation of the variational equation. 

Figure 1 

* 
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Since the shape of domain ft of the elastic solid is treated as the design 
variable, it is convenient to think of ft as a continuous medium and utilize 
the material derivative idea of continuum mechanics. The process of deforming 

ft to ft T by mapping ft^ = T(ft,x) may be viewed as a dynamic process of deforming 

a continuum, with t playing the role of time. A design velocity field can be 
considered as a perturbation of design variable (Refs. 2 and 4). 

Suppose z^(x ) is a solution of the variational equilibrium equation on 
the deformed domain ft T . Then the mapping z (x ) = z (x + xv(x) ) is defined on 
ft and z t (x t ) depends on x in two ways. First, it is T the solution of the 

boundary-value problem on ft T . Second, it is evaluated at a point x T that 

moves with x. Existence of the pointwise material derivative z is shown in 
Ref. 2. If z T has a regular extension to a neighborhood U T of the 
closure ft of ft T , then the partial derivative z' exists. One attractive 
feature of the partial derivative is that, with smoothness assumptions, it 
commutes with the derivative with respect to x^ (Ref. 2). (Fig. 2.) 


VARIATION OF DOMAIN 



z = — z T (x + xV(x) ) 


x=0 


= lim 
x-K) 


z^(x + xV(x)) - z(x) 


= z'(x) + Vz V(x) 



i = 1,2,3 


Figure 2 
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A common form of structural performance measure involves stress in an 
elastic solid. Consider a locally averaged stress functional ^ over a small 
subdomain flCA of the elastic solid, as shown in Figure 3, whePe g(a) is a 
stress measure such as von Mises stress or principal stress and m is a 
characteristic function that has a constant value on Q and its integral is 
1. The averaged stress measure depends on shape of thP domain in two ways; 
first directly on the domain over which the integral is carried out and second 
on the stress a that, in turn, depends on the displacement field z. 


Taking the first variation of ip , using material derivative formulas of 
Refs. 2 and 5, if*' is obtained. To obtain an explicit expression for in 
terras of the velocity field V, a variational adjoint equation is introfiuced by 
replacing z£ Z by a virtual displacement XGZ and equating terms involving X 
to the energy bilinear form a (X, X), yielding the variational adjoint 
equation for the adjoint variable X. 


STRESS SHAPE SENSITIVITY 


ffl n g(cr(z))da 
/// o g(°( z )) m da = 

H a da 

P 

3 • ■ . 

Iff Q [ I g i ,(z)o 1 J (z)]m da 

i,j=l a J p 

fff Q l [ l g i .(z)C 13 k£ (Vz k V )]m da 
i,j=l k,fc=l a 1J £ p 

+ ffl n g div Vm p da - /// fi gm p da /// fi m p div Vda 
3 

a 0 (*>*) = If [ I g . .(z)CT 13 (X)lm da, for all X £ Z 
a i,j = l a ij p 



Figure 3 
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Using the adjoint variable method of design sensitivity analysis (Refs. 2 
and 4) and the domain method of Ref. 5, an explicit and computable expression 
for ^ in terms of the velocity field V is obtained. Evaluation of the design 
sensitivity ty' requires the solution z of the original variational equation 

and the adjoint variable A of the variational adjoint equation. This is an 
efficient calculation, using finite-element analysis, if the original 
variational equation for z has already been solved, requiring only evaluation 
of the solution of the same set of finite-element equations with a different 
right side, called an adjoint load. 

For problems with smooth data in which stress is continuous, design 
sensitivity analysis results can be used for a pointwise stress functional. 

To obtain the formula, shrink the subdomain Q to a point x, where x£Q . In 
this case, the characteristic function become!? the Dirac delta measure.* 3 

Even though sensitivity analysis results for only a stress functional are 
presented here, the method is also applicable for displacement at a specified 
point x and eigenvalue design sensitivity analysis, as shown in Refs. 2 and 5. 
(Fig. 4.) 


r = *’(*) - a;( Z ,x> 

- /// l [ l g )]m dO 

n i , j =1 k , £= 1 a J F 

+ /// g div Vm dft - /// gm dft /// m div Vdft 


(z,X) = - Iff. J, [a lj (z)(vx i V.) + O^UXVz 1 V . ) ]dfi 
i,j = l J J 

3 . . 

+ fff 0 [ l a 1J (z)e J (X)] div Vdfi 

i.j = l 


3 T 3 

*'(X) = SSL l X^Vf 1 V)dfi + /// [ l f X 1 ] div Vd« 

V i=l 1=1 

+ //,{-! T i (VX lT V) + (7[ \ T^X 1 ] T n + H[ £ tV ] )(V T n) }df 


r i=i 


i=l 


1=1 


• Pointwise stress functional can be treated for problems with smooth data. 


Figure 4 
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For numerical implementation of shape design sensitivity analysis, the 
boundary r of the domain ft must be parameterized. There are several methods 
to parameterize the boundary T (Ref. 6). Since the result of shape 
optimization depends on the parameterization method used, it must be general 
and flexible enough to represent a large class of structural shapes. It is 
desirable that the parameterization method has the following properties: 
smoothness, fairness, required order of continuity, controllability in global 
and local senses, and a variation diminishing property. Among several 
parameterization methods, Bezier and B-spline surfaces are commonly used (Ref. 
6). Both Bezier and B-spline surfaces use a set of blending functions and 
are defined in terms of characteristic polyhedra. 

Points p (v,w), i = 1,2,3, on a Bezier surface are constructed by taking 

l 

linear combiniations of a set of blending functions B w (v) and D (w) and x, 

m,M n,N i 

coordinates c mnx of control points (vertices of the characteristic 

polyhedron). A iezier surface represented by a 4X4 array of points is shown 

in Figure 5. If a Bezier surface is used, positions c of the control 

mnx. 

points are shape design parameters. 1 

MODELING FOR SHAPE (BEZIER SURFACE) 



P v (v,w) 
X i 


M 


I 

m=0 


N 


I 

n=0 


c B (v)D (w), 
mnx^ ra ,M n ,N 


i = 1,2,3 


Positions of the control points are shape design parameters. 


Figure 5 
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The next step is to develop a general method of defining and computing a 
velocity field in the domain, in terms of the perturbations of the positions 
of control points. It is shown in Ref. 7 that regularity of the velocity 
field must be at least at the level of regularity oF the displacement field of 
the structure. This suggests use of displacement shape functions to 
systematically define the velocity field in the domain. Moreover, a velocity 
field that obeys the governing (elliptic) equation of the structure can he 
selected. That is, a perturbation of the boundary can be considered as a 
displacement at the boundary. With no additional external forces and a given 
displacement at the boundary, the finite-element equation can be used to Find 
the displacement (domain velocity) field, where {V^} is the given perturbation 
of nodes on the boundary, {Vd is the node velocity vector in the interior of 
the domain, and {f^} is the fictitious boundary force acting on the varying 
boundary. 

To use in Figure 4 for sensitivity computation, first perturb design 
parameter b i ^positions of control points), i=l, 2,***,k, a unit magnitude to 
obtain a boundary perturbation {V^}. Then domain velocity {V^} is obtained. 
Using {V^} and displacement shape functions, in Figure 4 can be evaluated, 
which gives 3^/ 3b. This method requires k solStions of the velocity 
equation. However, much as in adjoint analysis, this is an efficient 
calculation, requiring only evaluation of the solution of the same set of 
finite-element equations with a different right side for each unit 
perturbation of b^ , i=l,2, •••,k. (Fig. 6.) 

AUTOMATIC REGRIDDING FOR SHAPE DESIGN 


• Regularity of the velocity field must be the same as that of the 
displacement field 

• Use of displacement shape functions to define velocity field 

• Velocity field gives transformation mapping T(x,x) 



f K dJ 


• Solve above equation k-times 

• Excellent for boundary layer and/or substructuring technique 


Figure 6 
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The automatic regridding method presented in Figure 6 can be used with 
the boundary-layer approach (Ref. 8) and/or substructuring techniques very 
effectively. That is, if a large portion of the structure is fixed, except 
for the boundary layer (or substructure), then the part of {V^} that 
corresponds to the fixed portion can be set equal to zero, thus reducing the 
dimension of [K^]. 

Once a design change has been determined using an iterative design 
process, regridding of interior grid points can be carried out using {V^}. If 
the initial grid is optimized using an adaptive method (Ref. 9), the 
regridding method presented will tend to avoid distortion of the finite 
elements • 

To illustrate use of the automatic regridding method, a fillet problem 
(Figure 7) is used. In Figure 7, regridding is performed at three stages. It 
is interesting to observe that the method has a tendency to maintain 
orthogonality of the elements. That is, if the initial grid is regular, then 
the deformed grid tends to be regular. Also, the method presented can be 
utilized as mesh generator. That is, starting from a regular shape with a 
regularly patterned mesh (Figure 7(a)), the present method can be used to 
generate a mesh (Figure 7(d)) directly (Ref. 10). 

AUTOMATIC REGRIDDING FOR FILLET PROBLEM 






• The method can also be used as a mesh generator. 


Figure 7 
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To demonstrate use of the automatic regridding method for shape design, 
an engine bearing cap (Ref. 11), subject to oil film pressure and a bolt load, 
is treated (Figure 8). Oil film pressure is a radial pressure loading, 
assumed to be uniform. The engine bearing cap is modeled as a three 
dimensional elastic solid. Due to symmetry, only the right half of the cap is 
analyzed. The f inite-eleroent configuration and loading conditions are shown 
in Figure 8. The material used is steel, with Young’s modulus and Poisson’s 
ratio of E * 1.0 x 1(T psi and v = 0.3, respectively. The finite-element 
model shown in Figure 8 contains 82 elements, 768 nodal points, and 2111 
active degrees of freedom. For analysis and design velocity fields, the ANSYS 
finite-element STIF 95 (Ref. 12), which is a 20-node isoparametric element, is 
used. As in Ref. 13, implementation of design sensitivity analysis is 
performed outside the ANSYS finite element code. 

The shape design variables for this problem are: The shape of the 

varying surface T , distance Cc of clamping bolt center line AB, and distance 
c^ of edge from tne cap centerline. For surface T , a Bezier surface with a 
4x4 array of points is used. For simplicity, only X 2 ~coordinates of four 
control points c^ through c^ are allowed to vary. That is, surface T has 
curvature in the x^-direction only. 


ENGINE BEARING CAP 



• ANSYS STIF95 (20-Node Isoparametric element) 

• 82 elements, 768 nodes, and 2111 active DOF 


Figure 8 
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The expression for design sensitivity ij; f of averaged von Mises stress 

over individual finite elements is given in Figure 4, where g(cr) is von Mises 
1 ? 

stress* Define ^ and as the functional values for the initial design h 

P P 2 1 

and modified design b + 6b, respectively* Let At|> = and let be the 

P P P P 

predicted difference from sensitivity analysis. The ratio i^/Ai^ times 100 is 
used as a measure of accuracy; i.e., 100% means that the predicted change is 
exactly the same as the actual change. Notice this accuracy measure will not 
give meaningful information when A^ is very small compared to because the 
difference A^ may lose precision due Co the subtraction ^ 

Numerical result with a 1% uniform design change; i.e., 6b - 0.01 b, are 
shown in Figure 9 for randomly selected finite elements. Results given in 
Figure 9 show excellent agreement between predictions ^ and actual 
changes Aij^, except in elements ^ and 57. However, the magnitudes of actual 
change A^ are small for those elements. 


SHAPE DESIGN SENSITIVITY FOR ENGINE BEARING CAP, 6b = 0.01b 
(AVERAGED VON MISES STRESS OVER FINITE ELEMENTS) 


El. 

No. 

♦J 

< 

P 


(xp* / A' l> xl00)% 
P P 

1 

9829.4564 

9727.3229 

- 102.1335 

- 109.7298 

107.4 

5 

11444.4800 

11448.0190 

3.5390 

0.4482 

12.7 

10 

17933.5910 

17964.5170 

30.9260 

29.8750 

96.6 

14 

34270.5140 

34294.7650 

24.2510 

23.7614 

98.0 

20 

12670.2480 

12634.3500 

- 35.8980 

- 38.4216 

107.0 

26 

7311.4083 

6999.4094 

- 311.9989 

- 321.7022 

103.1 

30 

7234.2502 

7081.2085 

- 153.0417 

- 159.7947 

104.4 

35 

13328.4650 

13264.9790 

- 63.4860 

- 59.4243 

93.6 

39 

44231.0680 

42109.0220 

-2122.0460 

-2222.5504 

104.7 

44 

5998.6512 

5844.9335 

- 153.7177 

- 165.1199 

107.4 

48 

6822.9614 

6736.9477 

- 86.0137 

- 90.5011 

105.2 

53 

13634.1000 

12964.2560 

- 669.8440 

- 701.6882 

104.8 

57 

6121.4120 

6114.6667 

6.7453 

8.1242 

120.4 

62 

7041.7283 

6971.4204 

- 70.3079 

- 79.6051 

113.2 

66 

4787.5653 

4761.5085 

- 26.0568 

- 27.6278 

106.0 

71 

6541.8233 

6585.9308 

44.1075 

45.1422 

102.4 

75 

3820.6962 

3843.9362 

23.2400 

22.5210 

96.9 

80 

6240.3854 

6285.3485 

44.9631 

46.3209 

103.0 


Unit: psi 


Figure 9 
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A total hip reconstruction consists of a three-dimensional elastic solid 
composed of cement, a metal stem, and cortical and trabecular bone (Figure 
10). For simplicity, cortical and trabecular bone are modeled with the same 
material properties. Young's moduli! and Poisson's ratios for metal stem 
cement, and bone are: E 1 = 207 GPa, v 1 = 0.3, E 2 = 2.07 GPa, v 2 = 0.23, and 
E^ = 14.0 GPa, v3 = 0.3, respectively. 

The femur model shown in Figure 10 is obtained by approximating the real 
cadaver femur model of Ref. 14 with piecewise linear conical solids. For 
simplicity, structural and loading symmetries are assumed . Therefore, only 
half of the model is analyzed. A vertical load of 4000 N is applied at the 
tip of the metal stem. 

The finite-element model consists of 16 elements for the metal stem, 28 
elements for the cement, and 36 elements for the bone. ANSYS element STIF 95 
is used for all finite elements. The model has 525 nodes and 1335 active 
degrees of freedom. The model is assumed to be fixed at the distal end of the 
bone . 


TOTAL HIP RECONSTRUCTION (IMPLANT DESIGN) 




• Pointwise stress and strain energy density at interface. 

• 16 elements for stem, 28 elements for cement, and 36 elements for bone 
(all ANSYS STIF95) • 

• 525 nodes and 1335 active DOF. 


Figure 10 
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There are 16 shape design parameters: bj through bg are the radius of the 
metal stem and bg through bjg are the radius of the outer surface of the 
cement, at different locations along the center line. Thus, b i+8 - b.p 
i=l,2,...,8 is the thickness of the cement at those locations. The shape of 
the outer surface of the bone does not change. 

The principal stress is used as a design failure criteria for the metal 
stem and bone, whereas strain energy density is considered as the design 
failure criteria for cement. 

Shape design sensitivity results for pointwise principal stress in the 
stem at the stem-cement interface are given in Figure 11, for a 5% design 
change in design parameter b 5 . The pointwise stress is measured at a Gauss 

point (out of 9 Gauss points) on stem-cement interface of each stem finite 
element. 

Results presented in Figure 11 show excellent agreement between 
predictions tp and actual changes Ai|> , except in element 6. However, the 
magnitude of actual change A\|> is small compared to the magnitude of ^ for 
this element, so accuracy of the difference is questionable. P 

SHAPE DESIGN SENSITIVITY FOR IMPLANT DESIGN, «t> 5 = O.OSbj 
(POINTWISE PRINCIPAL STRESS IN THE STEM AT THE STEM-CEMENT INTERFACE) 


El. 

No. 


' |/ 2 Ai|» <|>* xl00)% 

P P P P P 


1 

' 65.75792800 

65.74896400 

- 0.00896400 

- 0.00875783 

97.70 

2 

77.13410600 

77.24745600 

0.11335000 

0.11608011 

102.41 

3 

58.03037400 

58.53323000 

0.50285600 

0.52206340 

103.82 

4 

77.00421000 

79.96762700 

2.96341700 

3.01203420 

101.64 

3 

151.71708000 

146.27679000 

- 5.44029000 

- 5.35753070 

98.48 

6 

234.54156000 

234.78980000 

0.24824000 

0.68237420 

274.88 

7 

288.65995000 

291.58509000 

2.92514000 

3.00576120 

102.76 

8 

149.94087000 

149.70614000 

- 0.23473000 

- 0.25492036 

108.60 

9 

20.76092900 

20.75818400 

- 0.00274500 

- 0.00277719 

101.17 

10 

6.23888850 

6.22105300 

- 0.01783550 

- 0.01811896 

101.59 

11 

3.99426970 

3.91787910 

- 0.07639060 

- 0.07985700 

104.54 

12 

6.25765390 

6.73601410 

0.47836020 

0.48739250 

101.89 

13 

15.90449700 

15.06538300 

- 0.83911400 

- 0.91444092 

108.98 

14 

23.77727200 

23.71259200 

- 0.06468000 

- 0.06987854 

108.04 


Unit : MPa 


Figure 11 
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Shape design sensitivity results for pointwise strain energy density of 
cement on the bone-cement interface are given in Figure 12, for a 5% design 
change in design parameter bg. The pointwise strain energy density is 
measured at one of the Gauss points at the bone-cement interface of each 
cement finite element. 

Results presented in Figure 12 show excellent agreeement between 
predictions i|»' and actual changes except in element 41. However, the 

magnitude of for this element is small compared to others. 

Even though results of sensitivity analysis of a pointwise principal 
stress in the stem and pointwise strain energy density in the cement are 
given, for variations of one design parameter for each, variations of all 
other design parameters yield similar results. Shape design sensitivity 
results for pointwise principal stress in the bone at the bone-cement 
interface and for pointwise strain energy density in the cement at the stem- 
cement interface are found to be excellent. 

SHAPE DESIGN SENSITIVITY FOR IMPLANT DESIGN, fibg = 0.05b 9 
(POINTWISE STRAIN ENERGY DENSITY IN THE CEMENT AT THE BONE-CEMENT INTERFACE) 


El. 

No. 



Ai|/ 

P 

% 

(f'/Aif/ xl00)% 
P P 

17 

2.693386 

2.864695 

0.171309 

0.185526 

108.30 

18 

1.324330 

1.346854 

0.022525 

0.025306 

112.35 

19 

1.358676 

1.373181 

0.014505 

0.016123 

111.15 

20 

2.968939 

2.965287 

-0.003652 

-0.003972 

108.76 

21 

6.532172 

6.527846 

-0.004325 

-0.004688 

108.39 

22 

6.197117 

6. 196119 

-0.000998 

-0.001068 

107.01 

23 

12.301795 

12.302323 

0.000528 

0.000569 

107.74 

38 

5.474089 

5.847445 

0.373356 

0.398447 

106.72 

39 

2.187812 

2.236401 

0.048590 

0.053682 

110.48 

40 

2.045186 

2.077065 

0.031879 

0.034058 

106.83 

41 

3.616023 

3.616629 

0.000606 

0.000478 

78.88 

42 

10.974028 

10.976652 

0.002624 

0.002725 

103.84 

43 

16.638003 

16.640659 

0.002656 

0.002837 

106.81 

44 

22.454411 

22.455967 

0.001556 

0.001666 

107.05 


Unit: kJ/m 3 


Figure 12 
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A doubly curvatured arch dam (Figure 13) that is similar to one treated 
by Wassermann (Ref, 15) is optimized using higher order finite-element 
approximation and the continuum shape design sensitivity analysis method 
presented here. 

The dam structure and loading conditions are assumed to be symmetric with 
respect to the crown cross section. Thus, only half of the dam is analyzed. 
Also, it is assumed that the dam foundation is rigid, and the gravel concrete 
is homogeneous and behaves elastically. Concrete’s elasticity modulus and 
Poisson’s ratio are E * 21.0 GPa and v = 0.2, respectively. Water and concrete 
weight densities are 10.0 kN/m 3 and 24.0 kN/m 3 , respectively. 

To parameterize two surfaces (water and free sides), Bezier surface 
parameterization is used with a 4x4 array of points. For a shape design 
parameter, the X 2 -coordinates of 32 control points are selected. The dam 
finite-element model contains 36 ANSYS STIF 95 elements, 315 nodal points, and 
726 active degrees of freedom. 

DOUBLY CURVATURED ARCH DAM 




36 elements (ANSYS STIF95), 315 nodes, and 726 active DOF. 


Figure 13 
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The principal stress is used as a design failure criteria. Principal 
stresses are measured at Gauss points on the surface of the dam (8 Gauss 
points for each finite element). Design sensitivity analysis results are 
tested for pointwise principal stresses. Excellent agreement between 
predictions and actual changes is obtained. 

The optimal design problem for the doubly curvatured arch dam is to 
minimize the volume of the dam, subject to constraints on pointwise principal 
stress on the surface of dam and thickness at the top of the dam. For 
iterative optimization, Pshenichny's linearization method (Ref. 16) is used. 
History of cost function and maximum constraint violation is shown in Figure 
14. After 17 design iterations, cost is reduced from an initial value of 
253,566 nr to 182,583 nr and the maximum tensile stress is reduced from an 
initial value of 3.084 MPa to 1.981 MPa. 

OPTIMIZATION OF DOUBLY CURVATURED ARCH DAM 


• Minimize volume subject to: 

Principal stress; - 10 MPa < < 2 MPa, i = 1,288 

Dam thickness; 6m < tj , j * 1,4 


ARCH DAM OPTIMIZATION COST-CONSTRAINT HISTORY 



Figure 14 
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A profile of the final design is shown in Figure 15. The final design 
shown in Figure 15 is rather different from Wasserman's design (Ref. 15), 
mainly in the bottom portion of the dam. The final design obtained here had 
developed a fillet in the bottom corner, which is not observed in Wasserman's 
design. 

In the crown cross section shown in the Figure 15, the middle portion is 
thinner than the top portion. From stress distribution in the final design, 
it is observed that the maximum tensile stress in this middle portion is well 
below the critical value of 2 MPa. Another interesting observation is that 
the compressive stress limit of -10 MPa has never been violated. In fact, at 
the final design, the maximum compressive stress is -5.202 MPa. 

A PROFILE OF THE FINAL DESIGN 




Figure 15 
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